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THE  MECHANISM  OF  THE  TRANSITION  FROM 
DEFLAGRATION  TO  DETONATION  IN  HIGH  EXPLOSIVES* 

Prepared  by: 

Carl  T.  Zovko 


Approved  by:  Evan  C.  Noonan,  Chief _ 

Physical  Chemistry  Division 


ABSTRACT:  Experimental  results  of  tne  study  of  spontaneous 

transition  from  deflagration  detonation  at  the  Naval 
Ordnance  Laboratory  Indicate  that  the  approach  to  the 
problem  can  be  In  two  stages;  the  first  Is  the  formation  of 
a  shock  from  pressure  waves  engendered  by  a  confined  deflag¬ 
ration,  and  the  second  the  3hock-inltlatlon  of  detonation. 
Since  a  preliminary  analytical  treatment  of  the  first  stage, 
reported  previously,  led  to  promising  results,  a  more 
extensive  IBM-704  program  has  now  been  undertaken.  Two 
numerical  codes  have  been  tested,  a  previously  developed  one 
based  on  the  so-called  "q-method"  and  a  special1  one  written 
for  this  program  which  avoids  amplitude  fluctuations  inherent 
In  the  “q-method"  and  thus  gives  a  more  realistic  representa¬ 
tion  of  a  shock  wave.  Representations  of  spontaneous  shock 
formation  obtained  uy  unc-  two  numerical  codes  and  by  the 
analytical  treatment  are  discussed  and  compared.  The  numeri¬ 
cal  methods  yield  the  temperature  as  a  function  of  time  and 
location  during  growth  of  tne  shock  and  thus  allow  a  study  of 
simple  cnemicai  kinetic  models.  Introduction  of  chemical 
kinetics  Into  tne  program  gives  a  basis  for  elucidation  of 
the  second  stage  of  the  transition  problem,  namely  shock- 
initiation  of  detonation. 
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Mechanisms  of  the  transition  from  deflagration  to  detonation 
in  explosives  and  propellants  are  presently  of  practical 
interest;  particularly  because  of  the  potential  hazard  posed 
by  large,  solid  propellant  rocket  grains.  Until  recently, 
little  was  known  about  the  processes  Involved.  Experimental 
studies  of  the  transition  were  published  in  NavOrds  5748, 
6104  and  6759.  This  report  is  concerned  with  theoretical 
interpretation  of  the  transition  process.  While  the  hydro- 
dynamic  model  is  one  dimensional,  considerable  insight  into 
the  mechanism  of  the  deflagration-detonation  transition  is 
provided.  It  is  part  of  a  continuing  investigation  in  this 
field . 

Tills  researcn  was  supported  by  NOL  Project  FR-59>  Transition 
from  Deflagration  to  Detonation. 

W.  D.  COLEMAN 
Captain,  USN 
Commander 

ALBERT  LJJGHTBODXy 
By  direction 
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THE  MECHANISM  OF  THE  TRANSITION  FROM 
DEFLAGRATION  TO  DETONATION  IN  HIGH  EXPLOSIVES 


INTRODUCTION 

While  the  phenomena  of  deflagration  (slow,  pressure- 
dependent  burning)  and  detonation  are  reasonably  well 
understood,  spontaneous  transition  from  one  regime  to  the 
other  is  still  in  early  exploratory  stages,  and  it  is  one 
of  the  major  unsolved  problems  in  explosives  technology. 
Gross  experimental  features  of  the  phenomenon  have  emerged 
only  recently  (1,2, 3, 4, 5).  It  appears  that  the  onset  of 
detonation  in  condensed  explosives  is  preceded  by  a  rela¬ 
tively  long  (up  to  80  psec  )  interval  of  rapid  burning  which 
propagates  at  a  fraction  (1/10  to  1/5)  of  the  steady  state 
detonation  velocity.  There  is  also  evidence  that  the 
actual  transition  from  rapid  burning  (sometimes  termed  "low 
order  detonation")  to  steady  state  detonation  takes  place 
rapidly  (within  several  microseconds)  at  a  plane  some 
distance  ahead  of  the  burning  front. 

The  evidence  thus  far  is  consistent  with  the  hypothesis 
that  the  onset  of  detonation  is  due  to  a  shock  wave  which 
arises  spontaneously  as  a  result  of  deflagration,  and  which 
initiates  detonation  in  unburnt  explosive.  The  hypothesis 
was  suojected  to  quantitative  scrutiny  at  this  Laboratory; 
in  addition  to  experiments  mentioned  above,  a  preliminary 
theoretical  treatment  was  carried  out  (6)  by  means  of  the 
following  model:* 

A  thermally  initiated  (slow)  laminar  flame 
progresses  into  a  homogeneous  solid  explosive 
charge.  Pressure  of  the  hot  products,  because 
of  rigid  confinement,  increases  steeply  and,  in 
consequence,  sends  compression  waves  into  unburnt 
explosive . 


*  Reference  6  gives  the  conceptual  and  analytical  basis  of 
the  computational  work  described  below,  and  it  will  be 
frequently  referred  to  in  the  subsequent  pages. 
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On  tnis  basis  it  was  shown  that  compression  waves  thus 
formed  coalesce  Into  a  shock  wave  within  10-15  cm  from  the 
region  of  thermal  initiation.  Since,  experimentally,  the 
typical  pre-detonation  distance  is  in  the  same  range  (6-14 
cm),  it  appears  reasonable  to  suppose  that  the  theoretically 
computed  snock  is  the  direct  cause  of  detonation. 

While  the  analytical  methods  thus  give  promising 
results,  it  is  very  desirable  to  extend  the  treatment  in  two 
ways:  first,  by  repeating  the  computation  using  different 
equation  of  state  parameters  and  different  shock-generating 
pressure  pulses;  and  second,  by  calculating  the  energy  (or 
temperature)  as  a  function  of  time  and  distance.  The  latter 
computation  can  then,  in  principle,  be  used  to  study  the 
chemical  kinetics  of  the  explosive  reaction  during  build-up 
and  thus  elucidate  the  transition  phenomenon.  Such  an  exten¬ 
sion  clearly  calls  for  machine  computation.  This  report 
gives  an  introduction  to  the  computational  program  which  is 
now  in  progress. 

The  report  consists  of  four  parts.  The  first  part 
describes  the  scope  of  the  program  treated  so  far  and  the 
equation  of  state  used.  The  second  and  third  parts  describe 
two  different  numerical  codes  for  the  IBM-704  computer  and 
compare  the  results  from  these  codes  with  the  previously 
obtained  analytic  results.  The  fourth  part  describes  the 
shock  formation,  chemical  kinetics  and  shock  initiation. 

SCOPE  OF  THE  PROBLEM 


A.  Hydrodynamics 

As  has  been  stated  above,  the  approach  to  the  problem 
of  transition  to  detonation  at  the  Naval  Ordnance  Laboratory 
has  been  via  two  stages.  The  first  one  is  formation  of  a 
shock  from  pressure  waves  engendered  by  a  confined  deflagra¬ 
tion.  The  second  one  is  shock-initiation  of  detonation. 

The  shock  formation  problem  is  programmed  in  the 
following  way:  The  difference  equations  for  conservation  of 
mass,  momentum  and  energy  are  written  down  as  applied  to  a 
one-dimensional  flow  problem.  The  explosive  charge,  which 
obeys  an  equation  of  state  described  below,  is  divided  Into 
N  zones  (0\N<500).  At  time  t  =  0  the  pressure  throughout  the 
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charge  is  fixed  at  a  low  but  finite  value  (P(t=0)  =  0.08  kbar). 
At  subsequent  times,  the  near  boundary  is  subjected  to  pre¬ 
scribed  pressures  increasing  with  time.  The  result  is  that 
compression  waves  of  increasing  amplitudes  travel  forward  from 
the  near  boundary. 

For  a  realistic  description  of  the  transition  process 
the  pressure  at  the  near  boundary  must  simulate  the  backing 
pressure  rise  in  a  confined  deflagration.  In  such  a  case  the 
theoretical  relationship  between  P  and  t,  derived  in  Ref.  6, 
is  given  by 

_  P 

t  =  K  dp  ,  (1) 

J  PQ  P(A-P)2 


where  K  and  A  are  constants;  at  low  pressures  this  is  ,  , 
sufficiently  well  approximated  by  the  exponential  P=pQe 
where  P0  (i.e.  pressure  at  t=0)  and  k  are  experimental  param¬ 
eters.  The  exponential  form,  which  was  used  previously  in 
the  analytical  treatment,  is  used  also  in  the  machine  compu¬ 
tations.  However,  an  indefinitely  long  exponential  pressure 
increase  would  be  unrealistic.  Decause  it  would  lead  to 
unreasonably  high  pressures  as  well  as  to  extremely  high 
values  of  dP/dt.  In  reality,  such  a  situation  does  not  occur; 
rather,  the  pressure  will  increase  until  the  confinement  is 
broken  and  then  decrease.  As  a  crude  simulation  of  such 
oehavior  tne  pressure  in  the  computation  is  allowed  to 
increase  exponentially  until  about  10  microseconds  after  the 
estimated  bursting  pressure  of  the  steel  casing  has  been 
attained;  thereafter,  the  pressure  is  assumed  constant.  The 
last  stipulation  may  be  at  least  partly  justified  if  one 
assumes  that  the  actual  pressure  decrease  is  relatively  slow; 
it  appears  rather  more  realistic  than  the  other  extreme, 
namely  a  discontinuous  pressure  drop  to  zero,  which  would 
cause  too  rapid  a  rear-rarefaction  to  set  in.  Thus  the 
assumed  near  boundary  condition  is 


P  =  P0ekt 

'  t  * 

^  =  pmax 

t  * 

^ (^max ) 
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Hydrodynamically ,  the  problem  of  coalescence  of  com¬ 
pression  waves  into  a  shock;  can  be  divided  in  two  parts. 

In  the  first  part  the  compression  is  isentropic  and  the  flow 
is  simple.  The  compression  energy  is  Es  and  the  temperature 
attained,  Ts ,  is  given  by  Ts  -  T0  =  Es,  where  Cv  is  heat 

capacity  of  the  explosive  and  the°c(mblent  temperature. 
This  part  was  treated  analytically  in  Ref.  6  by  the  method 
of  characteristics.  The  method,  in  fact  is  valid  only  for 
3uch  simple  flow  (i.e.  no  shocks);  it  does  not  give  a  basis 
for  further  calculation;  in  particular,  it  cannot  show  where 
and  when  the  shock  becomes  strong  enough  to  initiate  detona¬ 
tion  . 


The  second  part  of  the  problem  starts  with  the  overlap 
of  simple  waves.  The  flow  then  ceases  to  be  simple  and  there 
is  an  increase  of  entropy  across  the  compression  wave.  As  P 
and  dP/dt  at  the  near  boundary  Increase,  shock  compression 
conditions,  described  by  the  Rankine-Hugoniot  relations,  may 
ultimately  be  reached  at  the  front  of  the  disturbance.  Since 
the  most  important  part  of  the  problem  is  expected  to  be  the 
region  of  shock  formation,  i.e.  region  intermediate  between 
simple  flow  and  shock  conditions,  a  parameter  £  ,  whicn 

measures  the  extent  of  shock  nature  is  hereby  defined  such 
that 


El 


< 


(2) 


Eh"  Ej 

Here  E,  Ep  and  EH  are  actual  (computed),  isentropic  and  snock 

lion  energies  respectively,  corresponding 
(Since  simple  flow  is  isentropic 


(Hugonlot )  compress: 
to  a  given  pressure 

(  AS=0),  an  alternative  parameter,  0< 


AS  ^  1 ,  could  be 


defined  to  measure 
of  the  hydrodynamic 
£  >  0  respectively; 
corresponds  to  a  fu 


Sh 


the  extent  of  shock  nature).  The  two  parts 
problem  are  thus  characterized  by  £  =0  and 
the  upper  limit  of  the  parameter,  £  =  1, 

11  grown  shock. 


While  in  a  condensed  medium  the  difference,  for  a  given 
pressure,  in  energy  (and  consdquently  in  temperature)  between 
the  two  modes  of  compression  characterized  by  the  extreme 
values  of  £  will  r.ot  be  large,  the  difference  in  chemical 
reaction  rates  should  be  quite  considerable  and  may  mean  a 
difference  between  failure  and  initiation  of  detonation. 
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Hence  it  is  convenient  that,  in  addition  to  pressure, 
another  parameter  specifying  the  energy  be  known.  £  has  been 
chosen  because  it  gives  a  direct  indication  of  deviation  from 
simple  flow  conditions. 

B.  Equation  of  State 

A  generalized  Tait  equation  of  state  has  been  chosen 
to  represent  the  solid  explosive 

(P  +  B)  V  -  (P0  +  B)  VQ  =  (y-  1)  (E  -  E0)  (3) 

with  appropriate  values  of  the  constants  B,  and  y  ,  the  equa¬ 
tion  gives  a  remarkably  realistic  representation  of  the 
compression  of  solid  explosives  over  a  wide  range  of  pressures*. 
Combined  with  the  isentropic  condition, 

dE  =  -PdV,  (4) 

the  equation  reduces  to  the  form  used  in  Ref.  b  (which  does 
not  include  tne  energy): 


Explicit  equations  relating  the  various  properties  for  lsen- 
troplc  compression  and  for  shock  compression  on  the  oasis  of 
Eqn.  3  are  collected  in  Table  I. 

Tne  arbitrary  parameters  cnosen  in  Ref.  6  were  3  =  105 
kbar,  y  =  5.  The  choice  deserves  a  comment. 

If  Eqn.  (3)  is  to  oe  fitted  to  a  set  of  data  in  a 
certain  range  of  pressures,  tne  constants  B  and  yean,  in 
general,  be  assigned  any  convenient  values.  If,  however,  the 
lower  limit  of  the  range  is  P  =  o,  by  virtue  of  the  relation 

-  v  fW).  ■  l/TilT-B),  (6) 

*The  author  is  indebted  to  Dr.  S.J.  Jacobs  for  having 
pointed  out  the  promising  possibilities  of  this  extremely 
simple  equation. 
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the  value  of  B  Is  fixed  by 

Q  ^ 

B  =  22  (6*) 

•  o 

This  Is  certainly  the  case  in  the  shock  formation  problem, 
where  in  the  early  stages  of  shock  growth,  the  pressure  is 
quite  low.  The  value  B  *  105  kbar  used  in  Ref.  6  corresponds 
to  an  initial  sonic  velocity  Co  *  2.56  mm/psec,  which  is  an 
average  of  the  range  of  2.25  -  2.85  mm/psec  found  by  Majowicz 
(7)  for  a  series  of  explosives.  Thus  the  value  of  this 
parameter  is  realistic. 

There  is  no  doubt  tnat  the  value  of  7  =  5*  used  in  Ref.  6 
(and  by  some  earlier  workers),  is  too  low,  because  it  gives  an 
unrealistically  high  compressibility.  The  reason  why  the 
value  has  been  used  at  all  is  twofold.  First,  it  is  a  carry¬ 
over  from  calculations  of  high  pressure  gases,  such  as  detona¬ 
tion  products,  in  which  the  Eqn.  5  with  B  =  0  and  7  3  gives 

reasonable  results.  Second,  and  perhaps  more  important,  the 
choice  of  7  =  3  lends  convenient  tractabillty  to  hydrodynamic 
equations.  In  particular,  it  allows  the  boundary  path  in  the 
shock  formation  problem  to  be  evaluated  in  closed  form  (see 
Ref.  8);  this  would  be  impossible  for  any  value  7  > 3  (and 
probably  for  most  non-integral  values). 

Figure  1  shows  a  comparison  of  the  computed  P  -  V  relation 
for  two  different  sets  of  parameters  B  and  7  as  well  as  experi¬ 
mental  data  of  Majowicz  and  Jacobs  (9).  The  high  compressibil¬ 
ity  of  a  hypothetic  material  for  which  7  *  3  Is  evident.  The 
value  of  7  =  4.5,  on  the  other  hand  (combined  with  B  =  100  kbar) 
is  very  realistic,  and  it  is  the  current  choice  for  the  machine 
computations.  However,  since  the  analytical  treatment  exists 
(Ref.  o)  in  which  the  first  set  of  values  was  used  (7  =  3, 

B  =  105),  the  preliminary  computations  discussed  below,  were 
run  with  this  set  of  parameters  for  the  sake  of  comparison. 

The  sonic  velocity,  corresponding  to  B  =  100  kbar,  is 
2.5  mm/psec,  a  most  reasonable  value. 

NUMERICAL  SOLUTION  OF  HYDRODYNAMIC  PROBLEMS  (GENERAL) 

The  general  hydrodynamic  problem  is  a  solution  of  the 
equations  of  motion,  state  and  energy  release  subject  to 
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appropriate  boundary  conditions.  The  equations  of  motion  for 
a  one  dimensional  case  are 

po  Zli  =  _  ZZ  (conservation  of  momentum)  (7) 

dg  .  p  c>V  _  n  (conservation  of  energy)  (8) 

dt  ~  at  5t  “  u 


and  „n  dV  du  (conservation  of  mass),  (9) 

Ut  =  d* 

w lie  re  t  =  time 


p0  =  initial  density 
u  =  particle  velocity 
P  =  pressure 

E  =  specific  internal  energy 

^  =  heat  added  per  unit  mass  (from  chemical 

reaction) 

V  -  specific  volume 
X  =  distance 

x  =  Lagrange  coordinate  defined  by  the 

relation 


dX(x,t)  p0  V (x, t ) . 
d  x 

The  equation  of  state  is 


P  =  0  (E , V )  . 


(10) 


The  equation  of  chemical  energy  release  is 

dg  =  R(ti,E,V).  (11) 

at 
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One  way  to  obtain  a  solution  Is  by  numerical  techniques. 
This  consists  of  dividing  the  Lagrange  space  coordinate  (x) 
Into  a  number  of  equal  zones  and  approximating  the  differen¬ 
tials  in  equations  7,  8,  9  and  11  by  finite  difference  ratios. 
In  the  difference  equations,  the  dependent  variables  are 
usually  specified  at  the  interfaces  between  the  zones  or  at 
the  centers  of  the  zones. 

The  boundary  conditions  must  specify  the  values  of  the 
dependent  variables  for  all  values  of  the  Lagrangian  space 
coordinate  (x)  at  time  zero  and  for  the  end  points  (x  =  0  and 
x  =  xmax)  at  all  times.  Once  the  boundary  conditions  are 
specified  the  difference  equations  can  then  be  solved  to 
obtain  the  values  of  the  variables  at  the  interior  points. 

Most  differencing  scneines  have  the  limitation  that  they 
cannot  handle  discontinuties .  The  equations  of  motion 
(Eqns.  7,  8  and  9)  admit  discontinuous  solutions;  in  fact, 
the  discontinuities  are  the  most  interesting  parts  of  the 
solutions.  Two  methods  of  overcoming  this  difficulty  will  be 
discussed  in  the  next  section. 

The  time  increment  (4t)  used  cannot  be  chosen  arbitrar¬ 
ily.  A  stability  analysis  (Ref.  10)  of  the  problem  will  yield 
a  maximum  value  of  A t  with  which  reasonable  results  can  be 
obtained.  Stability  analysis  is  an  analysis  of  the  history  of 
an  arbitrarily  introduced  error.  Usually  a  critical  value  of 
A  t  will  be  determined  such  that  if  A  t  were  to  be  made  larger 
than  this  critical  value,  the  error  will  increase,  if  At  were 
to  be  made  smaller  than  this  critical  value,  the  error  will 
decrease  and  if  A t  is  made  equal  to  this  critical  value  the 
error  will  remain  constant. 

TWO  SPECIFIC  METHODS  OF  OBTAINING  NUMERICAL  SOLUTIONS  TO 
HYDRODYNAMIC  PROBLEMS 

Two  methods  of  handling  discontinuities  will  be  dis¬ 
cussed  in  this  section.  They  are  the  Richtmyer-von  Neumann 
"q"  method  (Ref.  11)  and  the  Lax  method  (Ref.  12).  In  both 
methods  discontinuities  are  approximated  by  steep  but  finite 
slopes . 

The  "q"  method  eliminates  discontinuities  by  the  inclu¬ 
sion  of  an  artificial  dissipative  term.  Physically  it  can  be 
considered  as  a  one-dimensional  viscosity.  This  dissipative 
term  "q"  is  defined  by  the  equation 
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_  (KAx) 
V 


2 


6u 

5x 


(12) 


and  (8)eb«ar:e''^'3peoUv3?sr  the  PreSSUre  (f)  and  **-•  <7) 


d  u 
d  t 


5  , 

“  qJ  (p  +  q)  and 


(15) 


%.  _  2S. 

at  at 


(p  ♦  q) 


av 


0 


(14) 


Equation  (9),  which  does  not  contain  pressure,  is  unchanged. 

/ 1  ,  \  Th®  resultant  set  of  equations  (Eqns.  (9).  Qo)  ] 

i  i  f?nd  do  not  ^ave  discontinuous  solutions  The  * 

solutions  of  the  modified  equations  and  h  or  S  L 

tions  are  very  nearly  the  same  except  in  region?  where  the 

Eq^ti?ns0(9fe(^)Si(n)eqy?^0nS/?S1  haVe  a  discontinuity. 
continuity  by ’asmooth^but  i^urie^  a™mate  «»  «•- 


The  Hugoniot  relation  across  a  shock. 


Ei 


-  E 


t  = 


-L 


(pi  +  pf)  (Vi  -  vf) 


(15) 


is  not  affected  by  the  inclusion  of  q. 


T^S  "q"  mefcilod  is  successful  in  that  it  eliminates 

aoiS2?iin?i  iS*  and  SiV?3  a  S°°d  appellation  to  the  true 
solution  in  every  aspect  except  details  of  the  shock. 

4  computer  (IBM  704)  program  which  utilizes  the  ”q" 

at  the  Naval  SS"1  "f'J  dl-«enenclng  scheme  was  constructed 

from  thlt  cod  L^oraCor?  °i'  w-  Walker.  Some  results 

this  code  will  be  discussed  in  the  next  section. 


considering  the  arbitrary  constant  as  the  product  of  K  and 
Ax  is  superfluous  at  this  stage  of  the  discussion.  How- 

Hiiffo  Then  the  dlfferential  equations  are  replaced  by  finite 
e<luatlous  the  A  x  mentioned  above  and  the/lx  used 
as  the  increment  of  the  Independent  variable  are 
K  is  a  dimensionless  constant  usually  nelr  unity.  1  31 ‘ 
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Another  method  for  handling  the  problem  of  disconti¬ 
nuities  was  devised  by  Lax  (Ref.  12).  While  the  "q"  method 
involves  a  quasi-physical  concept  and  a  modification  of  the 
equations  of  motion,  the  Lax  metnod  does  neither.  Rather, 
it  handles  discontinuities  by  the  nature  of  its  unusual 
differencing  scheme.  The  Lax  scheme  requires  that  all  of  the 
differential  equations  be  in  perfect  differential  form,  i.e. 


A 


dY 

5T 


dZ 

5T 


(16) 


where  A  is  a  constant.  This  differential  equation  1b  then 
differenced  in  the  following  way 


A 


,t+4t 


_  X  (yt  +  y 

2  '  x+4x 


1 

27T 


(Zx+4x  "  Zx-4x 


t  )) 

x-Ax'  / 

) 


U7) 


ddY 


In  so-called  normal  regions  where  ■  is  small, 

dx<= 


2  <Yx^x  *  Yx-<lx>  -  Yx  •  US) 

In  this  case  the  Lax  difference  scheme  approaches  an  ordinary 
forward  difference  scheme.  Therefore,  in  normal  regions,  the 
solution  obtained  by  the  Lax  scheme  approaches  the  analytic 
solution. 

d  2v 

In  regions  where  - 3  is  high  (at  shocks)  Eqn.  (18)  is 

dxd 

not  valid  and  the  Lax  scheme  comes  into  effect.  It  causes  any 
discontinuities  (or  other  extreme  changes)  to  be  replaced  by  a 
steep  but  smooth  change. 


As  mentioned  earlier,  the  equations  of  motion  must  be  in 
perfect  differential  form  if  the  Lax  scheme  is  to  be  used. 

The  equations  of  conservation  of  mass  (Eqn.  (9))  and  momentum 
(Eqn.  (7 )) and  the  equation  of  chemical  energy  release  (Eqn. 
(11))  are  already  in  perfect  differential  form.  A  fourth, 
independent,  perfect  differential  equation  must  be  constructed. 
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This  can  be  done  by  multiplying  Eqn.  (7)  by  u,  Eqn.  (8)  by 
p0  and  Eqn.  (9)  by  -P  and  adding  the  results 


Simplifying; 

p°  "ft  +  E  ^  35  “  -fe  (pu)  (2°) 

One  of  the  results  of  this  task  is  a  computer  (IBM  704) 
program  to  solve  hydrodynamic  problems  oy  the  Lax  method. 
Appendix  I  gives  a  description  of  this  program. 

comparison  of  Analytic,  "q"  and  Lax  Methods 

The  general  hydrodynamic  problem  solved  numerically  by 
the  "q"  and  Lax  methods  as  described  above,  will  now  be  com¬ 
pared  to  the  previously  obtained  analytic  solution  (6).  The 
same  equation  of  state  and  boundary  conditions  were  used  in 
all  three  calculations.  Equation  (2)  was  used  as  the 
equation  of  state.  It  was  assumed  that  no  reaction  took 
place  so  q(x,t)  was  set  equal  to  zero. 

The  following  boundary  conditions  were  used  in  all  three 
calculations 


p 

(o,t) 

=  P (0,0)  ekt 

for  t  ^  60  psec ; 

k 

=  .1 

psec-P 

P 

(0,t) 

=  P (0,60) 

for  t  60  |j.sec ; 

P 

^xmax* 

t)=  P(0,0) 

P 

(x,0 ) 

=  P(0,0) 

u 

(x,0) 

=  0 

P 

(0,0) 

=  .08  koars 

du. 

at 


-  u 


aE 

at 


QJL 

ax 


a 

at 


ay 

at 


-  ? 


-  p 


3.2 

a  x 


av 

a 
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V  (x,0)  =  .6245  cc/gra  and  E  (x,0)  =  2  x  104  ergs/gm  were  the 
values  calculated  1’or  an  adiabatic  compression  from  .001  to 
.08  kbars.  These  boundary  conditions  approximate  the  boundary 
conditions  realized  In  the  experimental  work. 

Figure  2  compares  the  P  -  X  plots  at  6b  psec  oDtalned 
from  the  analytic,  "q"  and  Lax  methods.  Except  for  fluctu¬ 
ations  in  the  plateau,  the  "q"  method  agrees  more  closely 
with  the  analytic  method  than  does  the  Lax  method. 

Figure  6  compares  P  -  X  plots  at  100  ^sec  obtained  from 
the  "q"  and  Lax  methods.  At  this  time  a  real  solution  would 
have  a  discontinuity  extending  slightly  below  the  plateau. 

The  "q"  method  gives  a  somewhat  closer  approximation  to  this 
discontinuity  than  does  the  Lax  method.  However,  the  "q" 
method  gives  severe  fluctuations  in  the  plateau,  while  the 
Lax  method  gives  none. 

Because  of  the  exponential  dependence  of  reaction  rates 
on  temperature  (l.e.  energy),  the  spurious  fluctuations 
inherent  in  the  "q"  method  render  the  "q"  method  almost  use¬ 
less  for  reaction  rate  studies.  Therefore  all  further 
numerical  work  discussed  in  this  report  is  based  on  the  Lax 
scheme . 

SHOCK  FORMATION  AND  INITIATION 

A.  Shock  Formation 

The  analytic  solution  (Ref.  6)  to  this  problem 
showed  that  a  shock  had  started  to  form  at  about  12  cm  from 
the  boundary  at  90  psec.  The  analytic  method  cannot  give  the 
rate  of  shock  growth. 

The  rate  of  growth  of  the  shock  is  illustrated  in  Fig.  4 
which  gives  energy-distance  (or  temperature-distance)  pro¬ 
files  at  several  different  times  as  computed  numerically  by 
the  Lax  method  using  the  previously  defined  boundary  condi¬ 
tions  and  equation  of  state.  The  generating  pressure  pulse 
was  allowed  to  Increase  exponentially  for  60  p3ec  so  that  the 
maximum  pressure  reached  was  Pmax  =  .5E.27  kbars;  thereafter 
the  boundary  pressure  remained  at  52.27  kbars.  In  Fig.  4, 
the  upper  horizontal  line  gives  the  energy  that  would  result 
from  a  shock  compression  to  .52.27  kbars;  the  lower  line  gives 
the  energy  that  would  result  from  an  isentroplc  compression  to 
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5 '*.2 7  kbars.  The  actually  computed  energy  has  increased 
perceptibly  above  the  limiting  isentropic  value  at  72.1  usee, 
at  which  time  the  compression  front  is  about  5.3  cm  from  the 
boundary.  However,  the  transition  from  the  isentropic  com¬ 
pression  to  the  shock  compression  is  continuous;  there  is  no 
sharp  point  of  shock  formation. 

The  growth  of  the  shock  is  also  shown  in  Fig.  5,  in 
which  the  parameter  £  (evaluated  at  the  compression  front)  is 
plotted  against  time.  The  figure  also  gives  the  location  of 
the  compression  front  as  a  function  of  time,  so  that  the 
extent  of  the  shock  nature  (  £  )  in  the  front  can  be  read  both 
as  a  function  of  time  and  distance. 

B.  Shock  Initiation 

Figures  4  and  5  show  that,  assuming  a  chemically 
inert  medium,  the  shock  wave  is  half  developed  (£  =  .5)  when 
the  compression  wave  has  travelled  about  16  cm  into  the  charge. 
The  next  step  was  to  see  whether  the  temperatures  generated 
were  sufficient  to  start  a  detonation  in  an  actual  (l.e. 
chemically  reactive)  explosive,  and  where  the  detonation  would 
start.  The  latter  point  is  of  particular  interest  because, 
experimentally  in  the  NOL-DDT  test  (Ref.  6),  the  detonation 
starts  about  15  cm.  into  the  charge. 

The  computed  point  of  initiation,  in  general,  could 
be  located  anywhere  between  the  boundary  and  the  wave  front. 

As  seen  in  Fig.  4,  the  layer  of  explosive  near  the  compression 
front  is  at  a  temperature  higher  than  that  of  the  boundary, 
but  its  residence  time  at  that  temperature  is  shorter.  The 
point  at  which  the  chemical  reaction  rate  becomes  sufficiently 
high  to  generate  a  detonation  wave  will  evidently  depend  on 
the  specific  parameters  used  in  the  computation. 

In  order  to  see  if  the  theoretical  model  agrees  with 
the  experiments,  a  simple  first  order  kinetic  model  was  used, 
i  ,e . 

Ba 

(l  -  F(x,t)J  Ae  (21) 


mass  fraction  of  burnt  explosive, 
preexponential  factor. 


c>F(x,t)  = 

at 

where , 

F (x, t )  = 
A 
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Ea  -  activation  energy, 

R  *  gas  constant. 


and 


T(x,t)  *»  temperature. 


T(x,t) 


is  defined  by 
T(x,t) 


the  equation 

-  t  + 

'  0  -v 


(22) 


Since  Q(x, t )  Is  the  energy  liberated  by  the  chemical  reaction 
at  a  point  (x,t)  and  H  Is  the  heat  of  explosion  of  the 
explosive,  then 


F  ( x ,  t )  = 

zlll 


Therefore,  equation  (21)  can 


o»Q(x,t) 

- 


be  rewritten  as 


4(x,t) 
A  H 


Ea 

HT (x, t ) 


(25) 


(«*) 


This  Is  the  explicit  form  of  equation  (11)  that  was  used  In 
the  following  calculations.  The  constants  In  equations  (22) 
and  (24)  are 

Cv  =  1.254  x  107  er8s/gmJc 

AK  =  5.016  x  1010  er88/gm 
A  =  10^  sec--*- 
Ea  =  55,000  ca^mole. 


A  computer  run  was  made  using  the  previously  discussed 
boundary  conditions  and  equation  (24)  to  compute  the  reaction 
rate.  The  maximum  pressure  was  52.27  kbars.  The  computed 
temperatures  were  too  low  to  cause  any  appreciable  reaction. 
The  run  gave  results  almost  Identical  to  the  run  represented 
In  Fig.  4. 

Figure  6  Illustrates  the  most  Important  result  of  this 
run.  It  Is  a  plot  of  F(x,t)  vs  x  at  several  different  times. 
It  shows  that  after  the  shock  is  partly  developed,  the  greater 
reaction  rate  In  the  Interior  (due  to  the  greater  temperature 
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increase  from  the  partly  developed  shock)  causes  the  reaction 
to  proceed  farther  than  it  does  at  the  boundary. 

A  subsequent  run  was  made  with  one  important  change. 

The  near  boundary  pressure  was, 

P(0,t)  =  .08  e1^  kbars  for  t  ^  67.5  psec 

k  =  .1  psec-^ 

P(0,t)  =  .08  e^*^  kbars  =  68.52  kbars  for 

t  ^67.5  usee 

The  higher  pressure  caused  the  temperature  (i.e.  energy)  to 
reach  higher  values  than  in  the  previous  run.  The  reaction 
rates  from  these  higher  temperatures  were  great  enough  to 
cause  the  reaction  to  go  to  completion.*  The  reaction  first 
went  to  completion*  17.7  cm  in  from  the  boundary.  Figure  7 
illustrates  the  course  of  the  reaction.  It  is  a  plot  of 
F  ( x ,  t )  vs  x  at  several  different  times. 

Figure  8  is  a  plot  of  T(x,t)  vs  x  at  two  different  times 
from  two  different  computer  runs.  In  one  run  the  material 
was  assumed  to  be  non  reactive;  in  the  other  run  the  material 
was  assumed  to  be  reactive.  The  times  chosen  were  slightly 
before  and  slightly  after  the  reaction  went  to  completion  in 
the  reactive  run.  The  interior  temperature  is  higher  than 
the  boundary  temperature;  this  coupled  with  the  exponential 
dependence  of  reaction  rate  on  temperature  caused  the  reaction 
to  go  to  completion  in  the  interior  before  it  went  to  comple¬ 
tion  at  the  boundary. 

Figure  9  is  a  plot  of  pressure  vs  x  at  several  different 
times  for  the  68.9k  kbar  maximum  pressure,  reactive  explosive 
calculation.  It  shows  the  development  of  the  detonation  wave. 


*  The  first  order  reaction  assumed  here  would  never  actually 
go  to  completion,  but  in  a  numerical  computation  the 
reaction  goes  to  completion.  The  error  is  completely 
negligible . 
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DISCUSSION 

Measurements  and  rough  calculations  (Ref.  6)  Indicate 
that  the  maximum  boundary  pressure  attained  In  the  NOL-DDT 
test  is  about  32  kbars.  Calculations  based  on  this  maximum 
pressure  and  homogeneous  first  order  kinetics  show  no 
appreciable  reaction.  This  is  not  surprising  since  it  is 
almost  certain  that  initiation  by  weak  stimuli  (e.g.  weak 
shocks)  requires  some  mechanism  of  stress  concentration 
(e.g.  occluded  grit  or  gas  bubbles).  The  important  result 
from  the  32  kbar  calculation  was  the  observation  that  the 
reaction  inside  the  charge  surpassed  the  reaction  at  the 
boundary. 

A  later  calculation  was  made  based  on  a  maximum  bound¬ 
ary  pressure  of  68  kbars.  This  was  done  to  compensate  for 
the  absence  of  stress  concentrating  mechanisms  in  the  model 
used.  This  pressure  was  adequate  to  cause  initiation  of  the 
explosive.  The  initiation  started  in  the  interior  at  a 
location  comparable  to  the  experimental  results. 

It  was  concluded  from  these  calculations  that  the  shock 
formation-shock  initiation  model  for  the  transition  from 
deflagration  to  detonation  is  sssentlally  correct. 
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APPENDIX  I 

As  a  part  of  this  task  a  computer  (IBM  704)  program  was 
constructed  to  3olve  hydrodynamic  problems  on  the  basis  of  the 
Lax  scheme.  The  program  consists  of  a  main  routine  in  which 
the  equations  of  motion  are  integrated  and  certain  other 
unchanging  operations  are  performed.  Calculations  involving 
the  equation  of  state  reaction  rates,  boundary  conditions  and 
stability  are  carried  out  in  subroutines.  Thus  if  any  of  these 
things  must  be  changed,  only  the  appropriate  subroutine  need  be 
reprogrammed . 

The  program  runs  according  to  the  flow  diagram  in 
Figure  10.  The  following  are  notes  to  Figure  10: 

(1)  The  stability  analysis  of  this  system  shows  that 
At  <  P0  A  x  — 

c 

t  is  computed  at  every  interface  and  the  smallest  value  is 
used . 


(2)  At  this  step  the  equations  of  conservation  of  mass 
and  momentum  are  integrated.  The  equation  of  conservation  of 
mass  (Eqn.  9),  when  differenced  according  to  the  Lax  scheme, 
becomes 


At 


(  \  (VX-MX  ♦  vx-4x> 


1 

Tlx 


,  ) 

x-Ax' 


Since  the  values  of  all  of  the  variables  are  known  at  t,  this 
equation  can  be  used  to  evaluate  Vt+  t.  Likewise,  the  equation 

of  conservation  of  momentum  (Eqn. 7*)  when  differenced  according 
to  the  Lax  scheme  becomes. 


Po 

IF 


^  yt+It 


(U 


x+Ax 


+  U 


t 

x-dx 


,) 


SIS  (P^+Ix  "  px-4x) 


Since  the  values  of  all  of  the  variables  are  known  at  t,  this 
equation  can  be  used  to  evaluate  U^+^t. 


Since  0  <  x  <  xmax,  the  above  two  equations  cannot  be 
used  to  evaluate  the  variables  at  t-Mt  where  x  -  0  or  xmax 
because  this  would  demand  values  of  the  variables  at  x  outside 
the  range  0  <  x  <  xmax. 
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(3)  At  this  step  the  equations  of  chemical  energy 
release  and  conservation  of  energy  are  Integrated.  The 
equation  of  chemical  energy  release,  when  differenced  accord¬ 
ing  to  the  Lax  scheme,  becomes 

If  \  ’  TZ  ((Wx  +  Qx-4x)  j  Rx 

Since  the  right  hand  side  of  this  equation  can  be  evaluated 

directly,  it  will  not  be  differenced  or  averaged.  r£  are 
evaluated  in  the  reaction  rate  subroutine.  Therefore,  this 

equation  can  be  used  to  evaluate  Qx+^* 

(4)  The  values  of  the  variables  at  the  boundaries  (i.e. 
at  x  =  0  and  x  *  xmax)  are  computed  in  the  boundary  value  sub¬ 
routine.  The  values  of  the  pressure  P  are  specified  by  the 
equations, 

P (0, t )  ■  .08  e0,lt  kbar  for  t  ^  tc0 

P(0,t)  =  .08  e0,11cokbar  for  t  t 

V*  u 

P(xmax,t)  *  .08  kbar  for  all  values  of  t 

The  values  of  the  other  variables  at  the  boundaries  are 
computed  from  P(0,t),  f(xmax,t)  and  equations  7,  9,  11  and  19. 

The  Lax  differencing  scheme  cannot  be  used  at  the  boundaries 

because  it  would  require  values  of  the  variables  outside  the 

range  of  0<x<x__  .  Therefore  a  different  differencing  scheme 
max 

is  used.  For  a  partial  differential  equation  of  the  following 
general  type. 


A 

A  W 

a 

(dZ) 

{dx}* 

the  following 

differencing 

scheme  is  used; 

A  / 

vt+dt 

yt 

Yt+dt  yt  \ 

m  [ 

I  - 

X 

I 

X 

x+/it  "  x+Ax  ) 

3  m  ( 

Z*  „  -  Zfc  +  Zt^t  -  Zt+At) 
x+Ax  x  x+Ax  x  / 

• 
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FIGURE  2 

Pressure-Distance  Profiles  as  Computed  by  the 
Lax,  "q"  and  Analytic  Methods 

The  time  of  the  computation  is  66.1  psec  after 
the  first  application  of  pressure. 
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FIGURE  5 

Pressure-Distance  Profiles  as  Computed  by  the 
Lax  and  "q"  Methods 

The  time  of  the  computation  is  100  p.sec  after 
the  first  application  of  pressure. 
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FIGURE  4  -  Energy  vs  Distance  at  Several  Different  Times  in  a  Non-Reactive  Material 


Irae 


FIGURE  6  -  Fraction  Reacted  vs  Distance  at  Several  Different  Times 


FIGURE  7  -  Fraction  Reacted  va  Distance  at  Several  Different  Times 


▼  a  i  a  j  v  .  x  v  u  u  » 


FIGURE  9  -  Pressure  vs  Distance  at  Several  Different  Times  with  a  Reactive  Material 
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